Mutations and intron polymorphisms in voltage-gated sodium channel genes of different geographic populations of Culex pipiens pallens/Culex pipiens quinquefasciatus in China

Background Culex pipiens pallens and Culex pipiens quinquefasciatus are the dominant species of Culex mosquitoes in China and important disease vectors. Long-term use of insecticides can cause mutations in the voltage-gated sodium channel (vgsc) gene of mosquitoes, but little is known about the current status and evolutionary origins of vgsc gene in different geographic populations. Therefore, this study aimed to determine the current status of vgsc genes in Cx. p. pallens and Cx. p. quinquefasciatus in China and to investigate the evolutionary inheritance of neighboring downstream introns of the vgsc gene to determine the impact of insecticides on long-term evolution. Methods Sampling was conducted from July to September 2021 in representative habitats of 22 provincial-level administrative divisions in China. Genomic DNA was extracted from 1308 mosquitoes, the IIS6 fragment of the vgsc gene on the nerve cell membrane was amplified using polymerase chain reaction, and the sequence was used to evaluate allele frequency and knockdown resistance (kdr) frequency. MEGA 11 was used to construct neighbor-joining (NJ) tree. PopART was used to build a TCS network. Results There were 6 alleles and 6 genotypes at the L1014 locus, which included the wild-type alleles TTA/L and CTA/L and the mutant alleles TTT/F, TTC/F, TCT/S and TCA/S. The geographic populations with a kdr frequency less than 20.00% were mainly concentrated in the regions north of 38° N, and the geographic populations with a kdr frequency greater than 80.00% were concentrated in the regions south of 30° N. kdr frequency increased with decreasing latitude. And within the same latitude, the frequency of kdr in large cities is relatively high. Mutations were correlated with the number of introns. The mutant allele TCA/S has only one intron, the mutant allele TTT/F has three introns, and the wild-type allele TTA/L has 17 introns. Conclusions Cx. p. pallens and Cx. p. quinquefasciatus have developed resistance to insecticides in most regions of China. The neighboring downstream introns of the vgsc gene gradually decreased to one intron with the mutation of the vgsc gene. Mutations may originate from multiple mutation events rather than from a single origin, and populations lacking mutations may be genetically isolated. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s40249-024-01197-1.


Background
Culex pipiens pallens and Culex pipiens quinquefasciatus are geographic subspecies of Culex pipiens, with Cx. p. pallens mainly distributed in northern China and Cx.p. quinquefasciatus dominating in southern China [1].The larvae of Cx. p. pallens often breed in polluted water bodies, and the adult mosquitoes are the main domestic nuisance mosquito species in China in terms of biting and blood-sucking, and this species is an important target for urban mosquito control.These mosquitos can also transmit bancroftian filariasis [2].
Insecticides are widely used to control mosquito vectors because of their high efficiency, rapidity, broad spectrum and convenience.Since the first discovery of resistance in 1908 [3], resistance has been found in an increasing number of insect populations.One of the more widely used insecticides, pyrethroids act on sodium channels; therefore, sodium channel mutations result in resistance to insecticides, that is kdr.Research has shown that a mutation at site 1014 of the SII6 fragment of the para-type sodium channel causes resistance in mosquitoes [4][5][6].Wild type L mutates to mutant type F [7]. Mutations at locus 1014 were found in both Cx.p. pallens and Cx.p. quinquefasciatus [7,8].The emergence of insecticide resistance in wild populations of mosquitoes is a complex process, which involves insecticide selection pressure during the evolutionary process from a wild homozygous (SS) to a wild/mutant heterozygote (RS), which then evolves into a mutant homozygous (RR) [9].Mosquito resistance is caused by not only insecticide selection pressure but also individual migration, genetic drift, and neutral evolution during the long-term evolution of mosquitoes, which may affect the emergence and development of resistance [9].Mutations at this locus induced by selective pressure from pyrethroid insecticides may lead to reduced nucleotide diversity in haplotypes [10,11].Variations in the neighboring intronic region of the gene can provide information about the long-term effects of pyrethroid insecticides on the evolution of mosquito populations [12].Genetic diversity analysis of this locus, including the intronic region, can aid in understanding population genetic characteristics and predicting the risk of the localized emergence of resistance.However, at present, there are large time intervals and regional limitations in the monitoring of resistance in Cx. p. pallens/Cx.p. quinquefasciatus.For example, a survey of Cx. p. quinquefasciatus resistance in some areas of Anhui Province was performed in the 1990s, but the survey of Cx. p. quinquefasciatus resistance in north-central Anhui Province was assessed only after a period of 20 years [13].To understand the current situation of Cx. p. pallens/Cx.p. quinquefasciatus in China, this study examined mutations in Cx. p. pallens/Cx.p. quinquefasciatus from 22 provincial-level administrative divisions.We analyzed the composition of downstream introns adjacent to the vgsc genes to understand the effects of pyrethroid insecticides on the long-term evolution of mosquito populations.This study provides insight into the impact of drug resistance on the long-term evolution of mosquitoes and provides a scientific basis for future directions of drug resistance monitoring.

Mosquito collection
Sampling was conducted from July to September 2021 in representative habitats in 22 provincial-level administrative divisions in China, including 28 populations (Table 1).The larval spoon method was used to collect larvae, and the mosquito trapping lamp method was used to collect adult mosquitoes [14].A total of 8000 mosquitoes were collected from the 28 populations.

Mosquito species identification
Morphological and molecular biology data were combined for mosquito species identification.Morphological identification was performed first.The ventral and inner leaves of the lateral plate of the male stem of Cx. p. pallens/Cx.p. quinquefasciatus are simple, with a wide leaf-like extension and a blunt end [2].Subsequently, cytochrome C oxidase I (COI) in the mitochondrial DNA was selected for molecular biology classification and identification.The 656 bp COI gene was amplified with the following primers: COI-F: 5'-GGC CAA CAA ATC ATA AAG ATA TTG G-3' , COI-R: 5'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3' [15].Each PCR mixture contained 12.5 μl of 2 × EasyTaqPCR SuperMix, 1 μl of each primer (10 μmol/L), 2 μl of template DNA, and ddH 2 O to a final volume of 25 μl.The thermocycler settings were as follows: 94 °C for 6 min; 35 cycles of 94 °C for 60 s, 51 °C for 60 s, and 72 °C for 60 s; and 72 °C for 8 min.

vgsc gene amplification
The DNA of a single adult Culex mosquito was extracted using a magnetic bead-based microtissue genomic DNA extraction kit (Bioteke Corporation, Wuxi, China).The whole-genome DNA was used as a template to amplify the structural domain IIS6 of the vgsc gene using the following primers for CD: CD1-F: 5'-GAC CTG CCA CGG TGG AAC T-3' , CD2-R: 5'-TTG GAC AAA AGC AAG GCT AAG-3' [9].The PCR mixture contained 2 μl of whole-genome DNA, 12.5 μl of 2 × EasyTaqPCR SuperMix (TransGen Biotech, Beijing, China), 1 μl of each of the 10 μmol/L forward and reverse primers, and ddH 2 O to a final volume of 25 μl.The PCR conditions were 94 ℃ for 3 min; 35 cycles of 94 ℃ for 30 s, 60 ℃ for 30 s, and 72 ℃ for 30 s; and 72 ℃ for 8 min.The samples were stored at 4 ℃.PCR amplification products were detected using 1% agarose gel electrophoresis, and samples with clear bands and no trailing were selected for forward sequencing.

Data analysis
Seqman [16] software was used to compare and analyze the peak maps of the vgsc gene, determine the allele and genotype of the test insect, and calculate kdr frequency.MUSCLE comparison [17] of sequences was performed using MEGA 11 [18].MEGA 11 was used to construct neighbor-joining (NJ) tree.GeneDoc [19] was used to plot sequence comparisons.The TCS network [20] was constructed using PopART [21], and the resulting where RR is resistant pure or resistant homozygote, and RS is a heterozygote genotype of wildtype and resistance genes.

Polymorphisms at the 1014 locus of the vgsc gene
A total of 1308 sequences of domain II of the vgsc gene were obtained in this study.Codon 1014 of structural domain II was analyzed, and six alleles and six genotypes were identified (Table 2).The first category of mutations was nonsynonymous mutations.Nonsynonymous mutations were categorized into two types: the mutation of leucine to phenylalanine (L1014F) and the mutation of leucine to serine (L1014S).The second category was the base double mutation of L1014F/L1014S, where leucine is mutated to both phenylalanine and serine at the same time (F1014S).The third category was the base synonymous mutation (L1014L).The six alleles obtained were TTA/L and CTA/L, encoding leucine; TTT/F and TTC/F, encoding phenylalanine; and TCT/S and TCA/S, encoding serine.The frequencies of TTT/F and TTA/L were greater, with frequencies of 47.17% and 43.58%, respectively.The genotypes included wild-type homozygous L/L, wild-type/mutant heterozygous F/L and S/L, mutant homozygous F/F and S/S, and mutant heterozygous F/S, with frequencies ranging from 37.77% (F/F) to 1.61% (S/L).

Spatial distribution of the genes
The TTT/F mutant allele was detected in all 28 populations (Table 3).The frequency of TTT/F ranged from 98.94% (YNJH) to 2.08% (JLTH).TCT/S was detected in the JXNC and SDLC populations at a frequency of 1.20%.CTA/L was found only in the SDLC population at a frequency of 2.38%.The F/S mutation was detected in 14 populations of ZJNB, HNSY, and SDZB, with frequencies ranging from 26.67% (ZJNB) to 1.19% (HNLY).The frequency of S/L was less than 10.00% in nine populations in which it was detected.In addition, F/F was detected in all 28 populations, with frequencies ranging from 97.87% (YNJH) to 2.08% (JLTH).Overall, TTT/F was the predominant mutant allele type (Table 4).
The allelic and genotypic compositions of geographic populations within the same latitudinal range were similar.kdr frequency increased with decreasing latitude.The geographic populations north of 38° N predominantly contained TTA/L and L/L, with the allelic mutation being TTT/F, and they had fewer homozygous mutations and a single mutation type.The frequency of mutant alleles increased in populations at 30°-38° N latitude compared to high-latitude populations, with the emergence of TCA/S and TCT/S mutant allele types and the mutant heterozygote F/S.Alleles and genotypes tended to be enriched in mid-latitude populations.The percentage of mutant alleles in some populations south of 30° N latitude was more than 90.00%, and these populations were dominated by TTT/F.Lower-latitude populations exhibited increased mutant purities, mainly F/F, and the alleles and genotypes were more homogeneous.Geographic populations with frequencies less than 20.00% were concentrated north of 38° N. Frequencies of 20.00%-80.00% were mostly found at 30°-38° N. Geographic populations south of 30° N latitude had frequencies greater than 80.00%.Within the same latitudinal range, the frequency was significantly greater in the large city populations than in the other populations.The frequencies of BJCP, HBWH, SHSJ, CQSPB, and GDGZ were 36.52%,93.24%, 86.84%, 90.48%, and 84.15%, respectively, which were greater than those of other populations at the same latitude (Table 5).

Characterization of neighboring downstream introns
Twenty-one haplotypes were obtained from 614 D2 sequences.Sequence analysis revealed that the mutant gene of Cx. p. pallens/Cx.p. quinquefasciatus was adjacent to a downstream intron in which the 5' end was GT, and the 3' end was AG (see Additional file 1).Overall, the total AT content was greater than the total GC content.The difference between AT and GC was most pronounced at the first base site, and A was predominant at all but the second base site (Table 6).These findings suggest that AT values increase in the downstream intron adjacent to the mutant gene of Cx. p. pallens/Cx.p. quinquefasciatus, and adenine (A) is preferred at the codon site.

Haplotype analysis of neighboring downstream introns to the mutant gene
The phylogenetic evolutionary tree of twenty-one intronic haplotypes showed that most of the intronic haplotypes were clustered in one large lineage, with H6 and H16 as independent branches.The intronic haplotype of TTA/L is widely distributed.Mutant allelelinked intronic haplotypes are distributed mainly within the same major lineage, in which the TTT/F intronic haplotypes H8 and H10 originate from the same branch and are sister groups to each other.The neighboring intronic haplogroups H1 and H5 of TTT/F and TTA/L are sister groups to each other.The neighboring intronic haplotypes H11 and H14 of TCA/S and TTA/L are sister groups to each other (Fig. 1).The haplotype network diagram shows that H17 is the dominant haplotype (Fig. 2).H17 may be the ancestral sequence from which the other haplotypes may have evolved.The mutant alleles TTT/F and TCA/S are derived from different TTA/Ls, with at least one sensitive haplotype on each branch.The mutant allele TCA/S has only one haplotype, H11, in the neighboring downstream intron, and TTT/F has three haplotypes, H1, H8, and H10, whereas the wild-type allele TTA/L has 17 haplotypes in the neighboring downstream intron (Fig. 2).It can be inferred that the neighboring downstream introns of the mutant allele are derived from the introns of the wild-type allele and are homologous to some of the wild-type introns.The introns of the mutant allele tend to stabilize, suggesting that the mutation is correlated with the introns.

Discussion
A total of three mutations were detected in this study: one was a nonsynonymous mutation, in which leucine L1014 was mutated to phenylalanine L1014F (TTA-TTT) or serine L1014S (TTA-TCA).The second was a base double mutation of L1014F/L1014S, in which leucine L1014 is mutated to both phenylalanine and serine F1014S.The third was a base synonymous mutation of L1014L (TTA-CTA).Consistent with the results of existing studies [22][23][24][25], the L1014F mutation is not conserved.In this study, the synonymous mutation CTA/L was detected only in the Liaocheng population in Shandong, China.Most of the mutations at locus 1014 of Aedes albopictus from the Anhui and Hunan Zhuzhou populations were CTA/L.It is possible that CTA/L is nonmutated in Aedes albopictus because locus 1014 has a CTA/L background.In contrast, the wild-type Cx.p. pallens/Cx.p. quinquefasciatus is TTA/L, and although CTA and TTA both encode leucine, CTA may be a mutation in Cx. p. pallens/Cx.p. quinquefasciatus.In the Zibo  population, the CTA and TTC allele genotypes were not detected, in contrast to the results of Wei's study [24].This difference may be due to the following reasons.The studies used different sampling habitats and sampling times.Convenient transportation, frequent transportation of goods, and well-developed tourism increase the chances of mosquitoes' passive movement, which increases the possibility of gene exchange between populations and changes the frequency of allele genotypes of the populations.And the populations had a low content of mutant alleles.
Resistance has developed in Cx. p. pallens/Cx.p. quinquefasciatus in most areas.The laboratory strain of Cx. p. pallens (BJ strain) had a resistance frequency of 15.96% and a mortality rate of 50.00% under exposure to a diagnostic dose of deltamethrin of 0.025%.The BJ strain was resistant to deltamethrin [26].Combining the BJ bioassay results and the frequency, it can be seen that field populations with a frequency < 20.00% have developed resistance to pyrethroid insecticides.The greater the frequency is, the lower the sensitivity to pyrethroid insecticides.This indicates that the light-colored Culex mosquitoes in most parts of the country have become resistant to the insecticide.Even populations with low frequencies may have developed resistance.When the frequency of vgsc alleles in Culex pipiens is low, cytochrome P450-mediated metabolic resistance dominates.With continued selection pressure, as the number of screening generations increases, the expression of P450 genes tends to stabilize, genetic differences between individuals in the population decrease, the actual heritability of drug resistance in the population increases, and vgsc mutations dominate.A combination of biological and molecular methods should be used to determine the level of resistance in populations, especially those with low frequencies, in a comprehensive manner.
Wild population resistance is regional in nature.The frequency increases as the latitude decreases.This study revealed that the mutation in Cx. p. pallens/Cx.p. quinquefasciatus in China exhibits significant spatial heterogeneity.The number and frequency of mutations showed significant geographic differences.Geographic populations with frequencies < 20.00% were mainly concentrated north of 38° N latitude.However, geographic populations south of 30° N latitude have frequencies > 80.00%.The populations south of 30° N latitude, such as the HNSY, YNJH, CQSPB and HBWH populations, had frequencies > 90.00%.No wild-type allele genotypes were detected in the HNSY population, with a mutation frequency of 100.00%.This could be because temperatures are higher at lower latitudes, and temperature is one of the most important factors affecting the growth cycle of mosquitoes and the level of breeding.Tire-causing Cx. p. quinquefasciatus were caught throughout the year in Hainan and Guangdong.In Jilin and Liaoning, the peak occurred in August, and after October, there were basically no catches of Cx. p. pallens [27].Mosquitoes overwinter differently, affecting the number of mosquitoes in the second year, which in turn affects the genetics of resistance.The more mosquitoes that overwinter, the more mosquitoes there will be in the following year.The greater the percentage of resistant mosquitoes that survive overwintering is, the greater the probability of resistance inheritance.Shandong's Cx. p. pallens overwinter at high temperatures and high humidity and in light-sheltered kilns, and the larvae cannot overwinter [28].Cx. p. quinquefasciatus in Jiangsu overwinters in underground garages, stairwells on the ground floor of buildings, and private houses on the outskirts of the city in different forms [29].On the other hand, Yunnan Province has a high incidence of dengue fever [30], and Chongqing had a dengue outbreak in 2019 [31].Therefore, the frequency, dose and duration of insecticide use are greater in these areas than in others.Insecticide residue levels are increased at mosquito breeding sites.The duration of mosquito exposure to insecticides increases, leading to increased mosquito resistance.The higher kdr frequency also converges across urban areas.The frequency of pesticide use is related not only to natural factors such as latitude and temperature but also to social factors such as economic, demographic.BJCP, SHSJ, and HBWH are large cities in terms of population, economy, and transportation.These cities are densely populated and serve as important transportation hubs with high logistics and people flow.Developed regions have a long history of insecticide use, and mosquitoes are exposed to insecticides for long periods of time, increasing cumulative toxicity.The frequency of insecticide use is greater in large cities within the same latitude.Mosquitoes are more likely to develop insecticide resistance.
This study only included mosquitoes collected in 2021 and lacked previous data.In recent years, due to its relationship with dengue fever, research on Aedes albopictus has been increasing, and research on Cx. p. pallens has correspondingly decreased.The main focus of this study was to investigate the haplotypes of adjacent introns downstream of locus 1014 using the appropriate sequences.The detected vgsc genes may not necessarily contain the complete introns required for research, making it impossible to obtain additional historical sequence data for Cx.p. pallens.Another approach is to compare functional genes (such as vgsc genes) with neutral markers (such as mitochondrial DNA), which may more accurately demonstrate the combined impact of demographic and selection factors on population structure.

Conclusions
This study described the spatial distribution of mutations in 28 populations of Cx. p. pallens/Cx.p. quinquefasciatus in China.The mosquitoes in most areas had developed insecticide resistance, and the degree of resistance was regional.The adjacent downstream intron of the gene was polymorphic and associated with the mutation.These results highlight the importance of monitoring the level of pyrethroid resistance in field mosquitoes, and limiting the spread of alleles will help to prevent the development of insecticide resistance.Establishing a long-term monitoring mechanism for insecticide resistance and understanding the development trend of insecticide resistance are helpful for implementing integrated mosquito control measures according to local conditions.

Fig. 1
Fig. 1 Haploid evolutionary tree of adjacent introns at position 1014 of the vgsc gene.Note: All unlabeled data are TTA/L.Pie charts indicate percentages at different latitudes.FJ182226.1 from NCBI, is the vgsc gene of Culex quinquefasciatus

Fig. 2
Fig. 2 Haploid network of adjacent introns at position 1014 of the vgsc gene.Note: The colored circles indicate the detected haplotypes.The black dots indicate deletions or haplotypes not actually observed in the sample.The circle size indicates the haplotype frequency.Consecutive vertical lines indicate the presence of mutations between haplotypes.The colored dashed lines indicate different clusters.Unlabeled lines are TTA/L

collection information Codes Location Sampling amount Acquisition time Collection method Latitude and longitude
haplotype network was plotted.kdr frequency was calculated as follows: kdr frequency = RR+ RS 2 N

Table 2
Detection results of target sites in the Cx.pipiens pallens/Cx.pipiens quinquefasciatus vgsc gene TTA, CTA, TTT, TTC, TCA, and TCT are codons.L, F, and S are concatenated bases representing leucine, phenylalanine, and serine, respectively.TTA and CTA encode leucine.TTT and TTC encode phenylalanine.TCA and TCT encode serine.n represents the total number of test mosquitos

Table 3
Allele distribution at locus 1014 in Cx. pipiens pallens/Cx.pipiens quinquefasciatus at different latitudes TTA, CTA, TTT, TTC, TCA, and TCT are codons.L, F, and S are concatenated bases representing leucine, phenylalanine, and serine, respectively.TTA and CTA encode leucine.TTT and TTC encode phenylalanine.TCA and TCT encode serine.n represents the total number of test mosquitos

Table 4
Genotype distribution at locus 1014 in Cx. pipiens pallens/Cx.pipiens quinquefasciatus at different latitudes p. quinquefasciatus varies with latitude, from single to diverse, before shifting to single.The genotypes are more monotypic north of 38° N latitude.The genotypes in the 30-38° N region are diverse.Populations south of L, F, and S are concatenated bases representing leucine, phenylalanine, and serine, respectively the mutant phenotypes of Cx. p. pallens/Cx.p. quinquefasciatus is increasing under the selective pressure of insecticides.Mutant genes can be stably inherited in populations, suggesting that genes may be conserved in the wild.

Table 5
kdr frequency of Cx. pipiens pallens/ Cx. pipiens quinquefasciarus in 28 geographical populations SS is a sensitive gene homozygous.RR 1 is a resistance gene homozygous (F/F and S/S).RR 2 is a resistance gene heterozygote (F/S).RS is a heterozygote for the sensitive gene and the resistance gene

Table 6
Sequence characteristics of downstream introns adjacent to the mutant gene T represents thymine.A represents adenine.G represents guanine.C represents cytosine.% Indicates the proportion of the nucleic acid in the fragment.a